Survival of the Fittest: Cox Regression
Capstone 6940
Intro to Cox Regression
Cox regression is a popular regression modeling method used to predict survival rates given certain covariates.
“Survival” can refer to the development of a symptom, time to relapse after remission, or as a time to death [1].
Cox regression model is based on the hazard function h_x(t) with covariates at time t given by:
h_x(t)=h_0(t)\exp(\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_n) h_x(t) = hazard \ function\ h_0(t) = Baseline \ hazard \ function\
\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_n = linear\ combination\ of \ covariates\ and\ their\ coefficients
Analysis: Corporate Survival
By Kayla Boyd
Data Source: Employee Turnover Data from Kaggle
| Variable | Description |
|---|---|
| stag | Experience (time) |
| event | Employee turnover |
| gender | Employee’s gender |
| age | Employee’s age (year) |
| industry | Employee’s industry |
| profession | Employee’s profession |
| traffic | From what pipelene employee came to the company. |
| coach | Presence of a coach (training) on probation |
| head_gender | Head (supervisor) gender |
| greywage | The salary does not seem to the tax authorities. |
| way | Employee’s way of transportation |
| extraversion | Extroversion score |
| independ | Independence score |
| selfcontrol | Self control score |
| anxiety | Anxiety score |
| novator | Innovator score |
{r, warning=FALSE, echo=F, message=FALSE}
library(readr)
turnover <- read_csv("turnover.csv")
View(turnover)
library(tidyverse)
library(dplyr)
library(corrplot)
library(RColorBrewer)
library(ggfortify)
library(riskRegression)
library(survival)
library(stringr)
library(zoo)
library(ranger)
library(ggplot2)
library(readxl)
library(MASS)
library(ADGofTest)
library(survminer)
library(car)
# Check for null and duplicate values
cat("Number of missing values :", sum(is.na(turnover)))
cat("Number of duplicates: ", turnover %>%
duplicated() %>%
sum())
# Remove duplicate values
turnover <- unique(turnover)
# Check distribution of data
turnover %>%
ggplot(aes(x = stag, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) + scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") + theme_minimal() +
theme(legend.position = "top")
# Check for censored observations
n <- dim(turnover)[1]
cat((n - sum(turnover$event))/n * 100, "% of observations are censored")
# Separate variable based on type
NUM_COLS <- c("stag", "age", "extraversion", "independ", "selfcontrol", "anxiety", "novator")
CAT_COLS <- c("gender", "industry", "profession", "traffic", "coach", "head_gender", "greywage", "way")
# Transform CAT_COLS into categorical type
for (COL in CAT_COLS){
turnover[COL] <- turnover[COL] %>% unlist() %>% factor()}Number of missing values : 0Number of duplicates: 13
49.82079 % of observations are censoredExploratory Data Analysis
Show Code
# Categorical covariates
# Selection of covariates (that are only discrete)
turnover.cat <- turnover %>%
select_if(is.factor) %>%
mutate(event = turnover$event)
turnover.cat %>% ggplot(aes(x = gender, color = factor(event),
fill = factor(event))) +
geom_bar(alpha = 0.5) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
Show Code
Show Code
Show Code
Show Code
Show Code
Show Code
Show Code
Show Code
# Continuous covariates
turnover %>% ggplot(aes(x = stag, color = factor(event),
fill = factor(event))) +
geom_histogram(aes(y = ..density..), alpha = 0.5) +
geom_density(alpha = 0.05) +
scale_color_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") +
theme_minimal() + theme(legend.position = "top")
Show Code
Show Code
Show Code
Show Code
Show Code
Show Code
Show Code
Show Code
# Other encoding of discrete variables
turnover.num <- turnover.cat[-9]
levels(turnover.num$gender) <- 1:length(levels(turnover.cat$gender))
levels(turnover.num$industry) <- 1:length(levels(turnover.cat$industry))
levels(turnover.num$profession) <- 1:length(levels(turnover.cat$profession))
levels(turnover.num$traffic) <- 1:length(levels(turnover.cat$traffic))
levels(turnover.num$coach) <- 1:length(levels(turnover.cat$coach))
levels(turnover.num$head_gender) <- 1:length(levels(turnover.cat$head_gender))
levels(turnover.num$greywage) <- 1:length(levels(turnover.cat$greywage))
levels(turnover.num$way) <- 1:length(levels(turnover.cat$way))
turnover.num <- as.data.frame(apply(turnover.num, 2, as.numeric))
turnover.num <- cbind(turnover.num, turnover %>%
select_if(is.numeric))
# Correlation matrix
corrplot(cor(turnover.num), col = brewer.pal(10, "BrBG"), method = "square", diag = FALSE)
Survival Plots
Now we will plot our categorical variables (gender, industry, profession, traffic, coach, head_gender, greywage, and way) using the Kaplan-Meier Estimation Curve. The Kaplan-Meier allows us to visualize how each category of survival probability differs from the others.
Show Code
Kaplan-Meier Curve Stratified by Treatment Groups
K-fold Cross Validation
Original c-index: 0.65 Mean cross-validated c-Index: 0.64 